# TODO: Add comment
# 
# Author: wuping
###############################################################################


dispersionIndex=function(x,y,indexName="poisson",sampleMethod="quadrat",alpha=0.05,lambda=NULL){
	x.mean=mean(x)
	x.var=var(x)
	n=length(x)
	pvalue=NULL
	lower=NULL
	if(sampleMethod == "quadrat"){
		if(indexName == "poisson"){
			I=x.var/x.mean
			#pvalue
			pvalue=pchisq(I*(n-1),df=n-1)
			if(pvalue>0.5)
				pvalue=1-pvalue
			attr(I,"pvalue")=pvalue
			#confidence interval
			upper=qchisq(1-alpha/2,df=n-1)/(n-1)
			lower=qchisq(alpha/2,df=n-1)/(n-1)
			attr(I,"confidence")=c(lower,upper)
			testName="Chi-squared test"
		}else if(indexName == "green"){
			I=(x.var/x.mean-1)/(sum(x)-1)
		}else if(indexName == "morisita"){
			I=sum((x-1)*x)/x.mean/(x.mean*n-1)
			#pvalue
			I.chi=(sum(x)-1)*I+n-sum(x)
			pvalue=pchisq(I.chi,df=n-1)
			if(pvalue >0.5)
				pvalue = 1- pvalue
			attr(I,"pvalue")=pvalue
			#confidence interval
			upper=(qchisq(alpha/2,df=n-1,lower.tail=FALSE)-n+sum(x))/(sum(x)-1)
			lower=(qchisq(alpha/2,df=n-1)-n+sum(x))/(sum(x)-1)
			attr(I,"confidence")=c(lower,upper)
			testName="Chi-squared test"
		}else if(indexName == "standardized_morisita"){
			I=sum((x-1)*x)/x.mean/(x.mean*n-1)
			#browser()
			#I=n*(sum(x^2)-sum(x))/(sum(x)^2-sum(x))
			Mu=(qchisq(alpha/2,df=n-1,lower.tail=FALSE)+n*(x.mean-1))/n/(x.mean-1)
			Mc=(qchisq(alpha/2,df=n-1)+n*(x.mean-1))/n/(x.mean-1)
			if(Mc>1 & I>=Mc)
				Ip=0.5+0.5*(I-Mc)/(n-Mc)
			if(Mc>I & I >=1)
				Ip=0.5*(I-1)/(Mc-1)
			if(1>I & I>Mu)
				Ip=-0.5*(I-1)/(Mc-1)
			if(1>Mu & Mu > I)
				Ip=-0.5+0.5*(I-Mu)/Mu
			I=Ip
		}
	}else if(sampleMethod == "distance"){
		if(indexName == "eToe"){
			re=1/sqrt(lambda)/2
			ra=mean(x)
			I=ra/re
			pvalue=pweibull(ra,shape=2,scale=1/sqrt(pi*lambda))
			if(pvalue >0.5)
				pvalue = 1- pvalue
			upper=ra/qweibull(alpha/2,shap=2,scale=1/sqrt(pi*lambda))
			lower=ra/qweibull(1-alpha/2,shap=2,scale=1/sqrt(pi*lambda))
			attr(I,"confidence")=c(lower,upper)
			testName="Weibull test"
		}else if(indexName == "pToe"){
			I=mean(pi*lambda*x^2)
			pvalue=pchisq(I*2*n,df=2*n)
			if(pvalue >0.5)
				pvalue = 1- pvalue
			lower=qchisq(alpha/2,df=2*n)/2/n
			upper=qchisq(1-alpha/2,df=2*n)/2/n
			attr(I,"confidence")=c(lower,upper)
			testName="Chi-square test"
		}else if (indexName == "hopkins"){
			#x should be event to event
			#y should be point to event
			I=sum(x^2)/sum(y^2)
			pvalue=pbeta(I/(I+1),n,n)
			if(pvalue >0.5)
				pvalue = 1- pvalue
			lower0=qbeta(alpha/2,n,n)
			upper0=qbeta(1-alpha/2,n,n)
			lower=lower0/(1-lower0)
			upper=upper0/(1-upper0)
			attr(I,"confidence")=c(lower,upper)
			testName="Beta test"
		}
	}
	cat(paste(sampleMethod,indexName,"dispersion index:",round(I,4),"\n"))
	if(!is.null(pvalue))
	    cat(paste(testName, "is used; p-value:", round(pvalue,6),"\n"))
	if(!is.null(lower))
	    cat(paste((1-alpha)*100, "% confidence interval (",round(lower,4),",",round(upper,4),") of random distribution \n"))
	return(I)
}
